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ABSTRACT 

Mass segregation stands as one of the most robust features of the dynamical evolution 
of self-gravitating star clusters. In this paper we formulate parametrised models of 
mass segregated star clusters in virial equilibrium. To this purpose we introduce mean 
inter-particle potentials for statistically described unsegregated systems and suggest 
a single-parameter generalisation of its form which gives a mass-segregated state. 
We describe an algorithm for construction of appropriate star cluster models. Their 
stability over several crossing-times is verified by following the evolution by means of 
direct A^-body integration. 

Key v^rords: stellar dynamics - methods: statistical - methods: A^-body simulations 



1 INTRODUCTION 

, Observations show quite often an increased concentration of 
' massive stars towards the centres of young star clusters (e.g. 
; ONC - Hillenbrand & Hartmann 1998; NGC 2157 - Fischer 
'. et al. 1998; NGC 3603 - Stolte et al. 2006). This tendency, 
' known as mass segregation, can be of different origin: Ini- 
, tial mass segregation is sometimes considered (e.g. Murray 
& Lin 1996, Bonnell & Bate 2006) as a consequence of the 
formation of massive stars preferably in the densest regions 
, (i.e. the cores) of the parent gas clouds. On the other hand, 
the process of mass segregation is also known to be one of the 
most robust features of the two-body relaxation driven evo- 
lution of self-gravitating star clusters (Chandrasekhar 1942, 
Spitzer 1969). 

Several approaches were developed to setup a star clus- 
ter in the state of mass segregation. Gunn & Griffin (1978), 
Capuzzo Dolcetta et. al (2005) and others based their setup 
on multi-component King models (King 1965, Da Costa 
& Freeman 1976) with stars separated into several mass 
classes which interact with each other via smoothed poten- 
tials. This approach relies on solving of non-linear set of 
Poisson equations, which is possible for limitted number of 
components. A multimass models of star cluster with exact 
energy equipartition in the core, which also leads to mass 
segregation, was introduced by Miocchi (2006). Another ap- 
proach used e.g. by McMillan & Vesperini (2007) relies on 
segregation produced by N-body integration of initially un- 
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segregated systems towards the segregated state, i.e. it is 
equivalent to a simple redefinition of time t = 0. 

In this paper we describe a new class of models of star 
clusters with continuous stellar mass distributions and a 
parametrised degree of mass segregation. The models are 
motivated by a study of the process of mass segregation dur- 
ing dynamical evolution of a self-gravitating cluster, which is 
brieffy described in the following Section. We show that mass 
segregation strongly manifests itself in the energy space. In 
Section O we introduce convenient characteristics of a sta- 
tistically described ensemble and derive their form for the 
unsegregated state. We further introduce in a heuristic man- 
ner an alternative, single-parameter form of these quantities 
that gives constraints on the distribution function of a mass 
segregated system. Afterwards, we describe an algorithm for 
construction of the corresponding star cluster. In Section [4] 
we demonstrate the stability of the models by means of A'^- 
body integrations. Finally, Section [S] contains our conclu- 
sions. 



2 MOTIVATION 

The standard scenario of the dynamical evolution of an 
isolated cluster is shown in Figure [1] The cluster in this 
example is initiated as an unsegregated Plummer model 
which is then integrated numerically with the NB0DY6 code 
(Aarseth 2003). We consider 20000 stars with masses in the 
range O.2M0 < m < 50Mq following a power-law mass func- 
tion with Salpeter index a = —2.35. The stars are treated 
as point-mass particles interacting solely by means of grav- 
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ity and, therefore, there is no intrinsic length-scale within 
the model. Hence, we introduce a characteristic length- and 
time- scale: 



ro = iGA^Vl-Etotl and to = rl'^/VCM: 



(1) 



by means of the cluster total mass, Afc, and the integral of 
the equations of motion, the total energy, -Etot- The results 
can be scaled to any length scale, provided the identities ([l]) 
between ro, to, Mc and -Etot are fulfilled. For a Plummer 
sphere, i.e. the initial state of the example model, the half- 
mass radius of the cluster is ~ 0.77ro and to corresponds 
to the crossing time. In the following we assume physically 
plausible stellar masses, although only ratios rm /Mc do mat- 
ter from the theoretical point of view. For definiteness, our 
'canonical' model presented in Fig. [1] has Mc — 132OOM0 , 
which for Vh = Ipc gives ro — 1.3pc and to = 0.2Myr. 

During the pre-core collapse phase of the cluster evo- 
lution, massive stars sink to the centre, forming a tightly 
bound core. This process is visible either in terms of the 
contraction of the inner Lagrange radii, or in terms of a 
decrease of the specific potential energy of massive stars. 
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is the potential energy of the i-th star. Compared to the 
latter quantity, the specific kinetic energy, K{muin) = 
E^VE'^ii of the same subset of stars shows a less pro- 
nounced evolution. At time t ~ 50to the cluster reaches the 
state of core collapse. From that point onward, strong few- 
body interactions between the core stars occur, leading to 
the formation of massive binaries carrying a considerable 
fraction of the cluster potential energy and, at the same 
time, to high velocity ejections of massive stars. This pro- 
cess stops further contraction of the Lagrange radii. How- 
ever, the potential energy of the subset of massive stars con- 
tinues to decrease. Both the potential and the kinetic energy 
show large variations due to the dynamical formation and 
destruction of binaries during the post-core collapse phase; 
the average kinetic energy of the massive stars starts to in- 
crease systematically due to the ejections. Note, however, 
that all stars are still kept in the computation. 

Another view of the redistribution of (potential) energy 
among the stars is presented in Fig. (2] where we plot the 
internal potential energy. 
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as a function of mass of an ordered subset of stars. 
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at four different times. In an initial, unsegregated, state, 
^^Hub should be proportional to the second power of A4*^l, 
which is clearly the case for the bottom line in Fig. [5] As 
time proceeds, the slope of the curve gets shallower, i.e. mas- 
sive stars hold an increasing fraction of the potential energy. 
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Figure 1. Evolution of characteristic quantities of an isolated 
cluster of 20000 stars. Top: Lagrange radii (0.5, 1, 2, 5, 10, 25, 
50, 75 and 90 per cent of Mc) are plotted with solid lines; dashed 
line indicates half-mass radius of a subset of massive (m > 5Mq) 
stars. Middle and bottom: the specific potential, U, and specific 
kinetic, K, energy of a subset of stars in terms of specific total 
energy, Etot = Etot /Mc. In all panels, solid lines represent quanti- 
ties related to the whole cluster (Afguj, = Mc), crosses correspond 
to the subset of stars with masses m > 5Mq {Mguh = 0.2Mc) and 
open squares represent characteristics of subset with m > 13Mq 
(A/sub = O.lAfc). Plots are obtained as an average over 100 runs. 
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Figure 2. Internal potential energy of a subset of stars as a 
function of its mass. From bottom to top the lines correspond to 
t = 0, 20, 40 and 60to- The plot represents an average over 100 
different realisations of the cluster presented in Fig. [l] 



At core collapse, the dependence can be approximated by 
another power-law function. 

The apparent monotonical evolution of the potential 
energy of stars of different masses during the whole course 
of the cluster evolution motivates us to parametrise the mass 
segregation in energy rather than in configuration space. 



This quantity implies that the mean potential energy of the 
i-th particle is 



(8) 



and the mean 'internal' potential energy of a subset of par- 
ticles is 



i-i 



(9) 
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The latter quantity will play an important role in the algo- 
rithm described below. In order to be unique, definition ([9} 
requires specification of the order of the particles in the set. 
Hence, we recall that we assume mi ^ m2 ^ ... ^ mjv. 

In order to take advantage of integral calculus, we will 
use replacements of summations; 



E 



nil. 
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Here, we use greek symbols to denote "continuous summa- 
tion index" . The equivalence in (|10p can be understood as an 
analogy to the discrete increment AM^^y^ = AA^, where 
A7V = 1. This trick introduces some error, in particular for 
small values of a and b and steep mass functions. Neverthe- 
less, it is useful to provide rather robust relations between 
individual quantities. 



3.1 Unsegregated state 



3 MODEL 

We consider an ensemble of N particles of masses rrii, i = 
1...N; we further denote with Mc = ^mi the total mass of 
the cluster. The state of the system is determined by spec- 
ifying N positions, r^, and conjugate momenta, pi, which 
altogether form a vector in a GA-dimensional phase space. 
For large N such a 'clean' state of the system is usually ei- 
ther not known or is not of particular interest. The system 
is then in a statistical sense conveniently characterised by 
means of a distribution function, _D]v(ri, pi, rjv, pjv), i.e. 
a probability density to find it in a particular state. For defi- 
niteness, we assume Dm to he normalised to unity. The mean 
value of an arbitrary physical quantity related to the system 
is obtained by integration over the whole phase space. 



(A) = J A{ri,pi, Tat, pjv)D]v(ri,pi, . 



, rjv,p]v)dr2 , (6) 



with dO, = d'^rid''pi...d''r]vd''p]v representing the phase 
space volume element. 

Specifying mean values of certain physical quantities is 
used to pose constraints on the form of the distribution func- 
tion in the case when it is not known explicitely. We assume 
the mean total energy (iJtot) is given. Restricting ourselves 
to systems in virial equilibrium, it follows that mean values 
of the total kinetic and potential energies are in balance, 
(A'tot) = -(-Etot), {[/tot) = 2{£'tot). We further assume 
that the system is characterised by the mean potential en- 
ergy between each two particles {i =^ j), 
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A commonly used scheme for construction of a cluster in a 
completely mixed (unsegregated) state is based on uncorre- 
lated drawing of positions, Vi, and velocities, Vi = pi/rrii, 
of individual stars according to a mass-independent single- 
particle distribution function /(r, v). This corresponds to 
the distribution function Dn in the form 

JV 

I>jv(ri, pi, rjv, pjv)dr2 = /(r,, Vi) d^r^d^v, (11) 

i = l 

with normalisation J /(r^, Vi) d^r^d^Vi = 1. For example, in 
case of a Plummer model. 
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for £ <0 and /(r, v) = otherwise. Here, Vp represents the 
characteristic radius of the Plummer sphere and 



GMc 



1 



(13) 



= 2"" rp ^1 + (r/rp)2 

is the specific energy of an individual particle. 

Regardless of the particular form of /(r, v), from sym- 
metries of the distribution function (|lip it directly comes 
out that the mean potential energy between two particles 
has to be a bilinear function of (and only of) their masses, 

(W^) = -Gm^m, f /('■Mv0.f(r.,v,) ^3^^^3^^^3^^^3^^ 



— —GGmiTrij 



(14) 

The integral in (|14p is an unknown constant C which is in- 
dependent of indices i and j. Its value can be easily obtained 
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by evaluation of (Ul^^^, defined in ((9)1, with the help of (fTO 



3.3 Building up the cluster 
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For i = N, i.e. M^^^^ = A/c, we require = (C/tot), 

which implies: 

2 {(/tot) 



C = 



GM2 



(16) 



For completeness, the mean potential energy of the i-th par- 
ticle, defined by formula (|8]) is 



^^^((7") = 2(t/tot) 



nii 



(17) 



3.2 Parametrisation of mass segregation 

From A'^-body models (Fig.[T} we see that it is predominantly 
the potential energy which is transferred between the light 
and massive stars, while their average kinetic energy remains 
nearly unchanged during the course of the cluster evolution. 
Hence, we will attempt to determine mass segregation in 
terms of mean potentials. In particular, we assume the mean 
inter-particle potential in the form: 

(t/'^) = 2{[/tot) ^^C/'^' (18) 

and consider several limitations to the term f/'-' : 

(i) it has to be symmetric with respect to the indices i 
and j. Only then will the total potential energy of the 
cluster be independent of the order of summation; 

(ii) it has to be positive and decreasing with increasing val- 
ues of indices i and j, so that massive stars (with lower 
indices) have lower specific potential energy; 

(iii) it should not depend explicitely on masses and rrij. 
Otherwise, core collapse could not be obtained for a 
cluster of equal mass stars. 



One of the simplest forms that fulfils these requirements is: 



2(1 -S)^ ([/tot) "^'"^^ /^-bM^, 



with S ^ being the tndex of mass segregation. In analogy 
to (|15|) . formula (|19|) implies 



and 



Ml 



(W) = 2(1 - S) {[/tot) 



Mc V Af< 



(20) 



(21) 



Clearly, 5 = corresponds to an unsegregated cluster while 
S > 1 would lead to a sign inconsistency of the potential 
energy of individual particles and (Utot)- Hence, only S £ 
(0, 1) should be considered as a reasonable value. The power- 
law form of (C/sub) (Mgub) is in accord with our motivation 
by the pre-core collapse evolutionary stages as depicted in 
Fig. [3 



Formula (|19p gives constraints on the distribution function 
of the cluster, although it does not determine it explicitely. 
The constraints expressed in terms of (U^) can, however, 
be used to construct0a corresponding star cluster by adding 
one by one the individual stars from the ordered set. 

The position of each added star is generated ran- 
domly (with isotropically distributed orientation) according 
to some 'underlying' distribution function n(r). The poten- 
tial energy of the (sub)cluster, U^, is calculated and com- 
pared with the desired mean value determined by eq. (fT9)l 
and @. (In the numerical code we drop the integral ap- 
proximation to calculate (U^) and evaluate it by means 
of summation in order to achieve better consistency.) If the 
difference | — {U^) \ is smaller than some given limi10, 
then we proceed to the next star in the set. Otherwise, we 
generate another position of the i-th star, until the match is 
adequate. 

The method for construction of the cluster described 
here has to be considered as a way how to find some state 
conforming to the given TV constraints. Hence, it is natural 
that the solutions do depend on the form of the underlying 
function used for generation of trial positions of added stars. 
The fewer trials are needed to find a matching position, the 
more likely is the final state close to a maximum of the 
distribution function Dn- By estimating contributions to the 
potential energy by individual particles (see Appendix E)) . 
we have found that a good underlying function (that needs 
on average less than 1.5 trials per particle) is given by: 



n(r) oc (rl{MU) + r^j 



-5/2 



with 

rp(ACb) = 



Ml 



GM^ 



\{Utot)\ 1-S \ Mc 



(22) 



(23) 



For S — 0, n(r) corresponds to the density of a Plummer 
model. Notice, however, that only for 5 = 0, the underly- 
ing distribution function is equivalent to the radial density 
profile of the cluster. The relation between the underlying 
distribution and the density profile of the obtained cluster 
is nontrivial due to the selection mechanism based on the 
check of Usuh vs. (f/gub) in each step. 



Velocity distribution 

At the end of the above procedure, positions and poten- 
tials of all particles are determined. In the next step we 
assign velocities to the stars such that the system is in a 
quasi-equilibrium state. As we assume the mean specific 
kinetic energy to be independent of the mass of the star, 
the distribution function of velocities has to be such that 
5 = — (Etot) /Mc. Furthermore, the velocity has to be 
corelated to the local gravitational potential, V{r), e.g. it 



^ A numerical C-code plumix for generating the cluster according 
to the algorithm described here can be downloaded from the Alf A 
web page: http;/ /www. astro. uni-bonn. dell 

2 In particular, for definiteness of the examples presented below, 
we considered | - (U^) | < | (U^) \l\Ji -I- 1 as a condition 
to accept the position of a particle. 
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has to fulfil V < ^2\V{r)\ in order to get a gravitationally 
bound system. 

We are motivated by the standard construction of the 
velocity distribution of a Plummer cluster (e.g. Aarseth, 
Henon & Wielen 1974), 



v{r)^q^2\V{r)\ 



(24) 



where q £ (0, 1) is a random number drawn from the distri- 
bution function 



(25) 



with mean square value, (g'^) — j.ln the case of our models, 
the explicit form of V(r) is not known, nevertheless, it can 
be replaced with /rrii which is calculated for each particle 
in the first step of the procedure. As for S 7^ 0, 



(26) 



is not independent of the particle index (mass), we cannot 
use directly the distribution (|25|) . Instead, we consider a 
generalised form 



n'{q;f3)^q^{l-q') 



2x/3 



(27) 



Then, (g'^) (/3) is a monotonically decreasing function, being 
singular at /? = —1. Hence, if we find /3 (see Appendix [B|) 
such that 



l-Etotl rrii 
(W) I Afc 



(28) 



formula (I24|) with q drawn from the distribution function 
(|27p will give the velocity of the i-th particle with mean 
square value 



2| jU') 
nii 



2\Et, 



(29) 



as required. For 5 = the method is equivalent to the stan- 
dard scheme used for the Plummer model. 
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Figure 3. Radial mean density profile of three different mod- 
els with S = (dotted), S = 0.25 (dashed) and S = 0.5 (soUd 
line). Clusters were built of 20000 particles with a Salpeter power- 
law mass function. In order to obtain smoothed density profiles 
even at very small radii, the profiles were obtained as an average 
of 20000 realisations with different initialisations of the random 
number generator. Half-mass radii of the clusters with different 
S are very similar (r^ 0.8ro; see Fig. [5] below) . Density is ex- 
pressed in units of po = Mc /rjj . 
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4 TESTS 

Fig. [3] shows radial density profiles of three models with dif- 
ferent values of the index 5*. In all cases, the clusters were 
built up with 20000 particles with masses according to the 
Salpeter power-law mass function used in Fig. [T] The case 

5 = corresponds to an unsegregated system. As the algo- 
rithm is very similar to a standard scheme for the Plummer 
model in this limit, the radial density profile obtains a char- 
acteristic shape with constant density core and outer parts 
with density falling as r~*/^. Increasing the index of mass 
segregation leads to considerable changes of the structure of 
the cluster. In the inner part (approximately up to the half- 
mass radius) the density can be approximated by a power- 
law, p(r) oc r-^-^^ and oc r'^ for S = 0.25 and S = 0.5, 
respectively. Note that in the latter case the density profile 
approximates the analytic solution of Lynden-Bell & Eggle- 
ton (1980) for core-collapsed (single-mass) star clusters (see 
also Baumgardt et al. 2003 for a numerical study of the pa- 
rameters of core collapse). 

Another view of the clusters' state, in terms of kinetic 
energy, is presented in Fig. [J] Here we plot the mean kinetic 
energy of stars within the inner region, r < 0.05ro, as a 
function of their mass. In the unsegregated state, S = 0, the 
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Figure 4. Mean kinetic energy in the core, rc < 0.05ro, as a 
function of the stellar mass. Models are identical to those in Fig.|3l 
i.e. dotted, dashed and solid lines correspond to 5 = 0, 0.25 and 
0.5, respectively. The kinetic energy was calculated for ten mass 
bins indicated with crosses. 

velocity distribution function is independent of the mass of 
stars and, therefore, K oa m everywhere within the cluster, 
including its core. For the sake of simplicity of the model we 
have posed a constraint of (K) oc m also for the mass seg- 
regated states. Nevertheless, according to eq. (|24p . the local 
mean kinetic energy depends on the particle mass. When 
placed at the same position, i.e. the same V{r), a light star 
will have on average a higher specific kinetic energy than 
a massive one as its velocity will be drawn from the dis- 
tribution (|27p with a lower /3, i.e. higher (3^)^ (note that 
index /3 of the velocity distribution depends on the index 
of the star, but it is independent of its position). Further- 
more, the selection criterion in the first step of the algorithm 
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Figure 5. Left panels: evolution of the Lagrange radii (like in Fig. [T] individual lines correspond to 0.5, 1, 2, 5, 10, 25, 70, 75 and 90 per 
cent of the cluster mass) of clusters generated with different values of the index of mass segregation. The right panels show the evolution 
of the half-mass radii of the whole cluster (solid), for a subset of stars with m > Mq (dashed) and a subset with m > 5Mq (dotted 
line). In order to distinguish trends from random fluctuations, the plots represent averages of 20 runs with identical values of the model 
parametres. 



allows low-mass stars to be placed in the innermost region 
extremely rareljjf], that is only if they hit the local minima of 
|y(r-)|. On the other hand, massive stars are allowed to en- 
ter local maxima of lV(r)| which stands as a multiplicative 
factor in eq. (|24|l . This effect slightly weakens the tendency 



^ For S = 0.5 stars from the lowest mass bin, m G 
(O.2M0, 0.35A/q), represent less than 0.03% of the total number 
of stars within O.OSro. 



towards energy equipartition in the core which, however, still 
remains a generic feature of the mass segregated models as 
it is demonstrated in Fig. |4l 

Dynamical evolution 

To test the stability of the models, we have used them as 
initial conditions for A'^-body integrations. Results are pre- 
sented in Fig. [5] The initially unsegregated system {S = 0) 
evolves rather smoothly without any apparent signs of in- 
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Figure 6. Mean kinetic energy in the core for the model with S = 
0.25. Dashed line represents the initial (already mass segregated) 
state, i.e. it is identical to the dashed line in Fig. |4] Solid line 
is a snapshot at t = 5to of that model integrated numerically. 
Dotted line corresponds to the model shown in Fig.[T]at t = 50to, 
evolving from an unsegregated state. Thin dotted line represents 
escape kinetic energy from the cluster centre. 



stability. Due to the high mass ratio, the process of mass 
segregation can be observed already on the time scale of a 
few crossing times. 

S = 0.25 gives a strongly segregated cluster. In this 
model the half-mass radius of stars heavier than SmM© con- 
tracts slightly, but it is already close to a saturated value. 
Note that its value is approximately one half of the half-mass 
radius of the whole cluster. This is in a good agreement with 
the model presented in Fig. [T] in the state of core collapse. 
In the right plot, which has a linear scale, we can see small 
initial oscillations of the half-mass radii which, however, are 
quickly damped. In general, this model can be considered as 
a quasi-stationary state very close to core collapse. 

The model with S = 0.5 shows more significant ini- 
tial oscillations as well as considerable overall expansion. It 
appears that for this value of 5", the algorithm produces a 
virially hot system with {Ktox) ~ 0.55(?7tot). This means 
that the criterion for accepting a newly added star on a 
particular position, which is formulated in terms of ((/^'ub), 
does not reproduce the mean potential energy (C/*) with 
sufficient accuracy. 

In order to test the stability of the models also in terms 
of the kinetic energy in the core, we show in Fig.|6]snapshots 
of the model with S — 0.25 at time t = and t = bto- 
We see that after a few crossing times the kinetic energy of 
the low mass stars settles at somewhat (approximately by a 
factor 1.3) higher values, while it remains nearly unchanged 
at the high-mass end. No further shift of the kinetic energy 
was observed for t > 5to. Interestingly, the new state, which 
settles after a few crossing times, fits very well to the state of 
a model which was followed from an initially unsegregated 
state to the state of core collapse (model form Fig. [1] at 
t = 50to). 

The escape kinetic energy from the core is linearly pro- 
portional to the stellar mass. Hence, in order to be bound 
to the cluster, light stars cannot have a kinetic energy equal 
to that of massive stars, i.e. the state of exact kinetic energy 
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Figure 7. Top: evolution of the Lagrange radii for an initially 
mass segregated (S = 0.25) cluster of 100000 stars with a Salpeter 
mass function. Bottom: half-mass radius of the whole cluster and 
subsets with m > Mq and m > 5Mq are plotted with solid, 
dashed and dotted lines, respectively. 

equipartition in the core is not possible even in the model 
with a rather high value of the index of mass segregation 
presented in Fig. (6] 

Mass function dependence 

The algorithm described in Sec. 13.31 does not depend ex- 
plicitely on the mass function (i.e. it can be used for any set 
mi). Nevertheless, its output is mass function dependent, 
as the constraints {U^^^,) depend on particular values of m;. 
We have performed tests with different sets in order to check 
the robustness of the algorithm. 

First, we considered the same mass function as before 
(i.e. m G (0.2, 50) and a = -2.35) but now with 100000 par- 
ticles. As we can see from Fig.[71 in terms of the characteris- 
tic radii the model with S — 0.25 evolves like its counterpart 
with 20000 particles, including the initial oscillations. A sim- 
ilar match was found also for different values of S, which we 
do not present here for the sake of brevity. 

A model with 20000 stars and a more artificial mass 
function (m G (O.IM©, lOM©) and a = —1.35) and mass 
segregation index S = 0.25 is presented in Fig. |8] Its be- 
haviour in terms of Lagrange or half-mass radii is similar to 
the case with the Salpeter mass function. The radial density 
and kinetic energy profiles are also similar to those presented 



8 L. Subr, P. Kroupa & H. Baumgardt 



10 



0.1 

0.01 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 




10 



4 6 



10 



Figure 8. Evolution of Lagrange radii (the same fractions as in 
the previous figures are considered) of a model with shallow mass 
function (o = —1.35, m £ (O.IA/q, IOMq)). Half mass radii in 
the bottom panel correspond to the whole cluster (solid), m > 
5M0 (Msub = 0.38Mc; dashed) and m > SM© (M^ub = 0.14Mc; 
dotted line) 



above. Hence, we can conclude that models with S < 0.5 are 
stable and their relaxational evolution is not influenced by 
apparently artifical effects. 

In both cases presented in this section, the few most 
massive stars played a less important role in the cluster 
dynamics, compared to the 'canonical' models presented in 
Figs. [3] - |6l These models helped us to reveal the origin of 
a slight flattening of the density profile, which can be ob- 
served in the case of a Salpeter mass function for S = 0.25 
(Fig. |3] dashed line). This effect is a demonstration of the 
dependence of (f/gub) on the masses of individual stars. In 
particular, {U^^^) determines the mean separation of the 
two heaviest particles, which is « O.OSro for a Salpeter mass 
function, but < O.Olro for the flatter one. Consequently, the 
mean density in this region stays approximately constant as 
the mass included is determined mainly by the two most 
massive stars. 



for construction of quasi-stationary models of mass segre- 
gated star clusters. For the sake of simplicity, we have per- 
formed tests with simple power-law mass function. Never- 
theless, the approach does not depend on a particular form 
of the mass function and the standard IMF (Kroupa 2001, 
2007) can be used as an input. Notice also that even for a 
cluster of equal mass stars the algorithm will lead to a sys- 
tem with a desired level of 'energy segregation' which is in 
general the process that drives the star clusters towards core 
collapse. 

Finally, let us remark that the index of mass segrega- 
tion is likely to be related to the entropy. For 5 = the 
system is highly symmetric in terms of (U^-''). On the other 
hand, in the limit of 5 = 1 it is required that all binding 
energy is carried by the two most massive particles, which 
is usually considered a state of maximal entropy of the self- 
gravitating system. We suggest that the statistical approach 
based on characterisation of the system by mean values of 
suitable physical quantities related to subsets of stars de- 
serves further investigation, providing us, hopefully, with a 
deeper understanding of the thermodynamics of star clus- 
ters. 
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APPENDIX A: UNDERLYING RADIAL 
DISTRIBUTION 



The mean potential energy of a Plummer cluster of mass Mc 
icteristic 

37r GM^ 



and a characteristic radius rp is 



32 



(Al) 



The contribution of mass from the interval {Msub,Afsub + 
dA/sub) is, approximately. 



d{t/sub) 



37r GM.ub dM.ub 
16 rn 



(A2) 



(this relation is exact only for rp = const). On the other 
hand, formula (1201) implies 

d ([/..b) = (2 - 25) ([/,o.) ( ' ^ . (A3) 



Afc 



Combining ^A2\ and (|A3|) gives an estimate (|23|l for 
rp(Afsub). 



5 CONCLUSIONS 

We have introduced a way of parametrisation of self- 
gravitating systems in terms of mean inter-particle poten- 
tials. We have demonstrated that this approach can be used 



APPENDIX B: PARAMETRISATION OF THE 
VELOCITY DISTRIBUTION 

The mean square value of a random number from an interval 
(0, 1) and probability density n' [q; (3) = (1 — q^)*^ is 
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(q^)^ = I^/Jp , (Bl) 
where 

Ip= f'q^l-qYdq (B2) 
Jo 

and 

Jp= r q^l^q^dq. (B3) 



Integrals (|B2|) and (|B3|) can be evaluated analytically for 
integer and half- integer /3 J5 — ^ by means of recursive for- 
mulae: 

/, = ^/.-r with 1-.,.-%, Io^\ (B4) 
and 

"^'^ ^ 5 + 2/3 "^""^ ^^^"^ "^"^''^ ^ I ' '^°^\' ^^^^ 

In order to find /3 giving (g^) according to equation (f28)) we 
start from /3 = —1/2 and evaluate recursively (9^)^ until 
upper and lower limits (g^)^j ^ {<f) ^ W^) 132 found. 
Then, we interpolate between /3i and (52- 

The Plummer model {S = 0) requires (g^)^ ~ \, 
i.e. /3 = I for all stars. Mass segregated models need 
{1^)0 < 3 fo'' massive stars and (9^)^ > \ for light 
ones. The procedure for finding appropriate (5 will fail for 
{q^) p > {l'^)_i/2 ~ I which, however, is not required even 
for the lightest star in the model with 5* = 0.5. 
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